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Abstract 


A moment method for computing 3-D radiative transport is applied to axisymmetric 
flows in thermochemical nonequilibrium. Such flows are representative of proposed aero- 
brake missions. The method uses the P-1 approximation to reduce the governing system 
of integro-differential equations to a coupled set of partial differential equations. A numer- 
ical solution method for these equations given actual variations of the radiation properties 
in thermochemical nonequilibrium blunt body flows is developed. Initial results from the 
method are shown and compared to tangent slab calculations. The agreement between the 
transport methods is found to be about 10 percent in the stagnation region, with the differ- 
ence increasing along the flank of the vehicle. 
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Nomenclature 


a Grid cell surface area (cm 2 ) 

B v Radiosity (W/cm 2 -s _1 ) 

C Constant in time term 

d Coefficients in Gauss-Seidel formulation 

G v Incident radiative intensity (W/cm 2 -s -1 ) 

I v Radiative intensity (W/cm 2 -s _1 -ster) 

H Emission coefficient (W/cm 3 -s _1 -ster) 

L Local error in Gauss-Seidel iteration 

n Surface unit normal 

q& Radiative flux (W/cm 2 ) 

r Relaxation parameter 

r Position vector 

R n Residual of Eq. 18 

s Path variable (cm) 

T Temperature (K) 

V Volume of grid cell (cm 3 ) 

x Cartesian coordinate (cm) 

y Cartesian coordinate (cm) 

z Cartesian coordinate (cm) 

, Computational variable for incident intensity 

e Surface emissivity 

( Third computational coordinate direction 

t] Second computational coordinate direction 

k' Absorption coefficient corrected for induced emission (1/cm) 

fj, Cosine of angle between 0 and z-axis 

v Frequency (1/s) 
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£ First computational coordinate direction 

r Optical variable 

to Solid angle (steradians) 

0 Direction vector 


Sub- and Superscripts 


i grid index in x or ^-direction 

/ cell center index in x or ^-direction 

j grid index in y or ^-direction 

J cell center index in y or ^-direction 

k grid index in z or ^-direction 

K cell center index in z or ^-direction 

m medium 

p Planck 

R radiative 

w wall 

v spectral 


n index of time step in iteration 


Introduction 

Figure 1 shows some important features of an aerobrake flowfield in a planetary atmo- 
sphere at the high speed typical of interplanetary flight. A significant portion of the heating 
experienced by such a blunt spacecraft can be due to radiation. The radiative heating can 
therefore be a strong driver on the design of the spacecraft heatshield. This heating is pre- 
dicted in two steps: first the radiative properties are computed at each point in the flowfield; 
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then the transport of radiation to the body surface is calculated. The first part of the prob- 
lem, for a flowheld including the features of Fig. 1, was treated in Ref. 1. For the radiative 
transport calculation the traditional model in a hypersonic, blunt body flow is the tangent 
slab or plane parallel approximation. 2-8 To further simplify the transport calculation, this 
model is sometimes applied to the visible and infrared region only, with absorption effects 
ignored. 9,10 

Recent work in aerothermodynamics 9 has provided some evidence that the ID tangent 
slab approximation may not be adequate even at the stagnation point for some of the vehicle 
shapes and applications under consideration for future space missions. It certainly will not 
suffice to compute radiative transport in the wake region where it is planned to locate 
unshielded payloads. No 3-D transport algorithm currently in use for flowheld solutions has 
been identified, though an effort toward this end is underway by Edwards et al. 11 . 

The objective of the present work is to initiate the development of a 3-D transport 
algorithm for radiative transfer in blunt body howhelds. This initial work is limited to 
axisymmetric forebody howhelds and does not include the effects of howheld coupling. The 
geometry is a deformed pie slice, with a cross-section like that shown in Fig. 2. Since an 
exact solution of the radiative transport equations in 3-D media is not feasible because it 
requires solving an excessive number of equations, an approximate method is used. The 
moment method is chosen here, reducing the problem from integro-differential equations to 
simpler partial differential equations. This method was originally developed for 1-D radiative 
transport, but recent work has extended it to 2- and 3-D geometries with various restrictions 
(see Ref. 12 and past work cited there). This paper discusses the application of the method 
to the current problem, including required changes for nonequilibrium, and presents some 
initial results. 


Radiative Transport Theory 

The complete equation of radiative transport has been derived in many standard texts. 
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Neglecting transients and assuming a non-scattering medium it becomes: 

dIA ‘' n) a) = ji(s) (i) 

Since for some frequencies k' v can be negative in part of a nonequilibrium flow, 1 it will 
pass through zero somewhere in the domain. The usual procedure of dividing by k' v to 
transform this equation to an optical variable would therefore introduce a singularity. For 
this work, then, Eq. 1 will be left in terms of the physical coordinate. 

Equation 1 applies at each frequency, z/, and for each direction, 0. A formal solution 
can be written using an integrating factor. The radiative transport can then be computed 
numerically for a finite number of frequency points and directions. Some method of select- 
ing appropriate frequencies at which to compute the radiative properties j e v and k' v in a 
nonequilibrium flow must be devised. That problem was discussed in Ref. 1. For the cases 
in this paper, 1066 optimized frequency points were used to compute the radiation occurring 

between 0.31 and 16.5 eV. 

Plane-Parallel Medium 


A common approximation for solving the radiative transport equation is to assume that 
the properties of the medium through which the radiation travels vary in only one direction. 
This situation is illustrated by the inset in Fig. 2. It may be a reasonable assumption in the 
stagnation region of a flowheld, where gradients are mainly in the normal direction. In this 
case s and 0 are replaced by z and fi (the cosine of the angle between 0 and the z-axis). 
The derivative d/ds then becomes fid / dz. To accommodate the boundary conditions in this 
case, / is split into forward (// >0) and backward (// <0) components. The formal solution 
for each can be written: 


I+(z, d) = 4 + (0, d) ex P - / <( 5 ", fi)/fids 

L Jo 

+ / jl{ s ' ■, d)/d ex P — / K 'A S ” ■> d)l dd s " 1 ds' 

Jo L Js' 


( 2 ) 


4 (a d) = h (a, d) ex p 



fi)/ fids" 
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fz r rz 1 

+ / h)/h exp - k' v (s" , fi) / fids" ds (3) 

•/^o L «/s / 

The z = 0 boundary is the vehicle surface, which is assumed to emit radiation according 
to Planck’s function at the wall temperature T w . The z = z 0 boundary is the freestream 
gas ahead of the vehicle which is assumed not to emit radiation. Along the surface normal 
(n = ±1), and introducing the optical variable: 

M z ) = f ^{s')ds' (4) 

Jo 

these then become: 

It{z, 1) = h P (T w )e ~ Mz) + [ Z me«W)dS (5) 

Jo 

K(z,-l) = r iKs'^-^ds' (6) 

J z 

To obtain the radiative flux in an absorbing, plane parallel medium it is easy to show 
that, if the exponential approximation is used to replace the exponential integrals, the above 
equations give the radiative flux, q v = + q~ } when t v is replaced by 2t v and a factor of 

7 r is added. The divergence of the radiative flux is obtained by differentiating this equation 
with respect to z, then integrating over the complete spectral range. The differentiation 
can be done analytically, or numerically from the solution for q v (z). The latter approach 
has been implemented here, following Nicolet. 2 The analytical approach requires a numerical 
differentiation of the integrand and would be of similar accuracy. 

Numerical Solution 

The integrals in Eqs. 4-6 require some care. Because k' v can be negative in some 
regions, 1 the log-linear variation assumed by Nicolet 2 in the RAD/EQUIL code cannot be 
used. Examination of k' v profiles suggests that a piecewise linear approximation is acceptable. 

For the intensity integrals, the flowheld grid spacing may be such that for strongly 
absorbing spectral regions A z is not small with respect to variations in the integrand. To 
integrate accurately, examination of t v profiles suggests a piecewise linear approximation be 
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used. Setting t v = cqs + A in each grid interval and approximating by its average at the 
two ends of the interval, the integral over the exponential term can be carried out. (The 
obvious approach of also setting = c 8 s + di has numerical problems for small 7y, resulting 
from terms which go to 0/0.) 

Three-Dimensional Medium 

When the tangent slab approximation is invalid, the directional variation of the radi- 
ation must be accounted for. The Modified Differential Approximation (MDA) of Modest 12 
is reported to give results of excellent accuracy for all optical conditions by splitting the 
radiative intensity into medium and wall components. The smoothly varying medium in- 
tensity can then be solved using the P-1 approximation. Adapting Modest’s derivation for 
the current non-gray, nonscattering, and nonequilibrium flow gives the following governing 
equation for the medium intensity (this result is equivalent to Eq. 1 when / is split): 

V • = jl - 0) (7) 

Taking the zeroth and first moments of this equation and integrating over 47 r steradians 
yields the applicable governing equations for the incident intensity, G v , and radiative flux, 
<fj,, of the medium: 


V -<Lvm = 47r 7Z - d v G Vm 

(8) 


(9) 


where the P-1 approximation has been used to simplify Eq. 9. 

To obtain the complete radiative transport solution, the contribution to q from the walls 
must also be considered. This is obtained from the following integral: 

( 10 ) 

7T J 47 r 

where the wall radiosity, B v , is given by the sum of emission at the wall temperature T w , and 
reflection of radiation emanating from other wall surfaces. In Modest’s analysis, the wall 
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reflection is assumed to be diffuse. The treatment of specularly reflecting wall surfaces adds 
considerably to the complexity of the problem and will not be addressed here. In any case, 
for a forebody flowfield the wall surface is convex so that no reflection of radiation emitted 
by the wall can occur (this simplification will not always apply in wake flows and so will be 
removed in future work on this method). Then the radiosity is just the emission from the 
wall: 

B v {r) = e u 7rI up [T w (r)\ (11) 

This contribution to the radiative flux must be computed for each grid cell in the medium 
by integration over all the wall surface elements to which it has a line of sight. 

If it is further assumed that the wall is cold, or that it makes a negligible contribution, 
then q w ss 0, and only the medium contribution need be considered. (Note that a cold 
black wall is assumed in many of the tangent slab transport analyses, since this gives the 
radiation incident on the wall. For generic studies in which the wall properties are undefined, 
this is a meaningful quantity for comparison.) This approximation is not unreasonable for 
aeroassist flowfields, since heating rates are relatively low in the upper atmosphere and the 
wall temperature will accordingly be relatively low during much of the flight. In fact, at the 
maximum wall temperature for Shuttle-type reusable thermal protection systems, the peak 
of the black body emission spectrum occurs at about 0.75 eV. For lower temperatures, the 
black body peak shifts to even lower energies. In this low energy spectral region the gas is 
quite transparent to radiation except in isolated atomic lines (see Fig. 3) so the wall can 
have little effect on the divergence of the radiative flux. In that case, V ' dR ~ V ' dm and is 
obtained by integrating Eq. 8 over frequency. Note that unlike the plane-parallel medium, 
numerical differentiation is not required here. 

Numerical Solution 

The Computational Fluid Dynamics (CFD) algorithm thus far used in this radiation 
study is a finite volume algorithm (the LAURA code of Gnoffo 13 ). To obtain a form of the 
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governing equations suitable for a finite- volume solution, Eq. 8 can be substituted in the 
divergence of Eq. 9 to obtain: 

VG^j = 3 k' v G v - 12 ? xjl 

In the hnite volume approach, this is now to be integrated over the volume of a single grid 
cell, V . By application of Gauss’ theorem (the divergence theorem), the volume integral on 
the left-hand-side can be transformed to a surface integral. Then 




J L \k 7 ^ Gi ' ' J da = Iv ~ 127r ^) ( 13 ) 

The integral on the left-hand-side is performed by assuming that the absorption coefficient 
k' v and the gradient of the incident intensity XjG v are constant on each cell face, leaving a 
summation over the six faces of the cell. On the right-hand-side, values at the cell center are 
used. The notation for the cell geometry is that of Gnoffo 14 , where /, J, and K denote cell 
centers, and z, j, and k denote cell faces. Then 


+ 

+ 


-T\7G„ 


* + l 


^i-\~ 1 . \/Giy I * Cli 

\ K , . 


J J,K 


K , \/G v j • dj+i f , \/G u j • a, 

v ' j+l \ v ' 3 J i t K 

"V^l ' 1 ( r yffi/ ] ■ Ufc 

< Jk+ 1 \< J k J/,j 

= V {^G V G V - 12? t] p v ) i jk 


(14) 


The remainder of the development will be restricted to axisymmetric flows, in which the 
second term in square brackets is zero. Cylindrical coordinates are not used, since the 
eventual goal is treatment of 3D flows. 

For the nongray gas found in a shock layer, the radiative properties k' v and j e v vary 
over orders of magnitude at various locations and frequencies. To minimize the numerical 
difficulties that this variation can introduce, some kind of normalization is desired. Define 

7 „ = ( 15 ) 

^cells 
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and 


— _ cells \ K v 
V ~ ^cells 

Both these quantities are constant for each frequency considered, and provide a measure of 
the approximate magnitude of the emission and absorption at that frequency. Figure 3 is an 
example of the variation of these average coefficients and demonstrates the wide variation 
in radiation properties for different radiative frequencies in a single flowheld. To use these 
normalizing quantities, Eq. 14 is divided by j v on both sides, and terms containing G v are 
multiplied by 1 in the form of Tc^/Tci,. A new variable , v = G u /(j v ~Ku) is defined. Eq. 14 then 
becomes: 



+ 


V, , ' a *+i ( / V? v I • t 

- K v ) 8 + 1 \ K v ) i J K 


Kv ^ \ _ / K v - \ 

! \7 5 V I ' a k+ 1 I ! \7 5 V j a k 

K v ) k+1 \ K v ) k 


= V 3k v k , v , v - 127ri^ 


1 v / I : K 


(17) 


This formulation reduces the variation of the unknown and coefficients of the equation rel- 
ative to Eq. 14, making the numerical solution easier. To further normalize the coefficients, 
this equation is multiplied on both sides by k' . 

The flowheld grid is not orthogonal, so the expression of the gradient in finite differ- 
ences requires a transformation to an orthogonal computational domain as shown in Fig. 4. 
Transforming to the new (£, () domain, and expanding the gradient of , v in second order 
central differences allows Eq. 17 to be expressed as: 


^1’ £V+1,A'+1 A ^2? vi+i,k A ^3? I 

+ ^4, Vi t K+l A ^ k • Vj t K A ^6? 1 / / /; _ 

+ ^7, A ^8’ v I — l,K A ^9 5 vi — i,K — l 

= —127 tV ik J -^- (18) 

where the coefficients are combinations of the geometry and the radiative properties 
resulting from collecting terms after Eq. 17 is differenced. 
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Equation 18 is a matrix equation for the unknown , Vl K of the general form Ax = b. 
The numerical solution of this type of equation can be obtained by a number of standard 
methods. The line Gauss- Seidel algorithm has been selected and applied along each normal 
grid line to capture the dominant gradients in the radiative properties. This algorithm 
results in a tridiagonal matrix equation to be solved for each normal grid line. The solution 
of such an equation can be obtained using the Thomas algorithm 15 with specification of 
appropriate boundary conditions as discussed below. To ensure convergence for this problem 
underrelaxation is necessary. It is introduced by defining a corrected value for the update of 
, v at the n + 1 time level as 

n+l' _ n , / n + 1 _ n I n+1 

’ vi t K ’ V I,K ‘ V i '/,if ’ v II( ) W )i i/ II( ~ 5 i/ II( 

where r is the relaxation parameter and is less than one for underrelaxation. The selection 
of r is discussed below. , ” +1 is substituted from Eq. 18 to obtain a new form of the matrix 
equation for , . Additional details of the derivation and the form of the coefficients 

are given in Ref. 16. 

So far this development has ignored the singularity introduced in Eq. 12 with the division 
by k' v . Two problems arise in this connection. The first is a loss in diagonal dominance in 
Eq. 18 resulting from the change in sign of k' v across a cell face. The second is the singularity 
itself. 

To obtain convergence at these frequencies, a time variation term, C dG v / dt, is intro- 
duced on the left-hand-side of Eq. 12. At the converged steady state this term is zero, 
but during the iterative solution it adds an additional term to the diagonal coefficient: 
d 5 => d 5 + CTy. The term CTy, also appears on the right-hand-side. The variable (7, 
which is related to a time step, can then be chosen such that diagonal dominance is recovered. 

To avoid the singularity in Eq. 12, note from the original equation, Eq. 9, that when 
k' v = 0, XjG v = 0. To enforce this in the numerical solution, k' v = 0 is assumed to occur at 
the cell face. The appropriate terms from Eq. 17 are then dropped for cell faces where k' v 
crosses zero. 
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This approach allows solutions to be obtained for all frequencies with a single algorithm: 
C is set to zero when k' v is everywhere positive, and to a value ensuring diagonal dominance 
when it is not; and the appropriate terms in the coefficients are set to zero where k' v 
crosses zero. The choice of the value of C is discussed below. 

Boundary Conditions 

Boundary conditions in the 3-D medium are illustrated in Fig. 2. The treatment of 
each is discussed below. 

Reflection boundary conditions should apply at the symmetry boundary. This boundary 
is a singularity in the computational grid, however, where one surface of a grid cell has zero 
area. Referring back to Eq. 17, this means the second term on the first line disappears. The 
remaining ^-derivative terms at the / grid line are replaced by first order forward differences 
to avoid differencing across the singularity. 17 The appropriate form of the coefficients is 
given in Ref. 16. 

For the medium only intensity, the wall boundary behaves as a cold wall. The Marshak 
cold wall boundary condition for the P-1 method is: 

2 [— - ll q v ■ n + G v = 0 (20) 

.Cjj J 

Replacing q v in this equation using Eq. 9 leads to a differential equation in the single unknown 
G v , which can be transformed to an equation for , v . This equation can be differenced at 
the wall boundary by recognizing that the grid lines are normal to the body surface. The 
V, v ■ n term is then simply the gradient of , v along the grid lines. This difference is best 
expressed in physical coordinates with a first order forward scheme, and can be incorporated 
in the Gauss-Seidel solution scheme with underrelaxation. The details are given in Ref. 16. 

Along the outflow “boundary”, the value of , v must be approximated. This is an 
arbitrary computational boundary through which radiating gas and radiation (inwardly and 
outwardly directed) both pass. The flowheld boundary condition used in LAURA is a zeroth 
order extrapolation (i.e., the derivatives of the flow variables are assumed to be zero across 
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this boundary). This boundary condition is known to provide better stability for CFD 
computations than other possible formulations, and so might be adopted for the radiation 
method. Several alternatives are possible at this boundary, such as a constant slope of 
, V (I,K) or constant radiation intensity. These various boundary conditions are found to 
have only a local effect on the solution. The selection of the “correct” boundary condition 
is hampered by grid distortion and a lack of resolution in the shoulder area of currently 
available CFD solutions. The choice is therefore made on the basis of stability rather than 
accuracy. The boundary condition that is applied in this work is 

(^(outflow boundary) = 0 (21) 

Though incorrect, this condition enhances stability while providing the correct trend of a 
fall-off in G v around the shoulder. Future work with complete flowhelds will move the outflow 
boundary downstream in the wake, where its influence on the flowheld will be much smaller. 

The freestream “boundary” is transparent and only outgoing radiation is considered. 
The freestream is assumed to be non-emitting. With these assumptions, this boundary can 
be modeled as a cold wall with complete absorption = 1.0 in Eq. 20). Differencing using a 
central difference for the {(-direction and a first-order backward difference for the {(-direction 
allows this equation to be written in the line Gauss-Seidel form. Underrelaxation can then 
be added as well. The final expression is reported in Ref. 16. 

For the freestream boundary at the axis of symmetry, a reflection boundary condition 
for G l/I _ 1 F is used. Using a first order forward difference at this point was found to lead to 
nonphysical results. 

Radiative Heat Flux to Wall 

The radiative flux to the wall due to emission in the medium can be found from solving 
Eq. 9 for ([„. The value of G v is obtained from the solution of the finite volume problem, , j,, 
multiplied by the normalizing factor (j ^v)- To obtain the flux directed to the wall, the dot 
product of q v with the wall-directed surface normal is required. Again recalling that the grid 
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lines are normal to the body surface, the gradient of G v can be expressed with a one-sided 
difference along the normal grid lines. As above, this is best done in physical coordinates: 


1_ (dG v \ 

^ V dTl J 

The total heat flux to the wall at each grid location is then obtained by a numerical inte- 
gration of Eq. 22 over the spectral range considered. 



Convergence Criterion 


The numerical solution of this set of difference equations requires the selection of an 
appropriate criterion to test for convergence. Because the radiative properties vary over 
many orders of magnitude at a single frequency, the usual error criteria must be modified. 
A local error function is defined as: 


_ \R n \ 
mean(, 


( 23 ) 


where R n is the residual of the difference equation, Eq. 18. The mean of , ” is the average 
of the minimum and maximum values of , ^ R for all / and K on the radiation grid. (For 
frequencies where k' v changes sign, it is defined as half the maximum value of , This 

definition reduces the error in regions where , ^ r is small, relative to the error that would 
be computed from the residual alone. Though the solution may not be completely converged 
at locations where , ^ r is small, these locations contribute little to the coupling with the 
flowheld. Therefore, the lack of complete convergence at these locations is deemed acceptable 
to reduce the computation time. 

Each frequency is converged individually, since the rate of convergence depends on the 
magnitude of the optical depth and varies widely for different frequencies. A typical conver- 
gence history for a single frequency is shown in Fig. 5. For other frequencies convergence 
requires 50 to several hundred iterations, and for a few frequencies even thousands of itera- 
tions. 
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Selection of Solution Parameters: r, C 


The relaxation parameter required to obtain a stable solution varies from frequency to 
frequency. No simple dependence of r on any radiative property has so far been discovered. 
The value r=0.5 has therefore been selected as a compromise, since it is found to provide 
generally good stability and convergence rates for all frequencies. 

The constant C in the time term added for frequencies with k' v crossing zero must be 
selected to recover diagonal dominance. Since C is related to the inverse of the time step, 
it must be kept as small as possible or the solution does not progress. To minimize the 
additional logic involved, an initial guess is made which will set d 5 = 1.1 d 5 . C is then 
successively doubled if tests indicate lack of diagonal dominance. 

It should be noted that the objective of this study is to demonstrate the feasibility of 
the method; not to optimize the numerical solution. Further improvements to this logic are 
almost surely possible and will be pursued as the method is used. 

Results and Discussion 

The LORAN radiation method 1 is used to obtain all the radiation predictions shown in 
this paper. Infrared, visible, and ultraviolet radiation from 0.31 to 16.5 eV is computed on an 
optimized spectrum of 1066 points. 1 LAURA flowheld solutions provide the nonequilibrium 
gas conditions. Coupling effects are not included in this study. 

Mars Return 

A flowheld solution has been obtained for one of a number of possible flight conditions 
identified for the return from a mission to Mars. 18 It consists of a 60° sphere cone with a 
1.08 m nose radius hying at 80 km altitude with a velocity of 12 km/sec. Results have been 
obtained for this case using both tangent slab transport and the MDA method. 

Figure 6 compares the wall radiation hux predictions from the two transport methods. 
The MDA result is found to be consistently lower than the tangent slab result. This results 
from a combination of geometric effects and howheld property gradients. Figure 7 presents 
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the ratio of the MDA to the tangent slab solution. It shows that the tangent slab solution 
applies best at the stagnation point (within 6 percent of MDA), and becomes progressively 
worse along the flank. The odd behavior at the shoulder is a result of grid problems in 
the flowheld solution used to obtain these results (these problems have since been fixed 
in LAURA). The qualitative trends of the two methods are similar, showing the radiative 
flux increasing with distance from the sphere-cone juncture. This increase is unlike the 
behavior of convective heating on such a body and results because radiation is a volumetric 
quantity. As suggested in Fig. 4, the standoff distance is much larger on the flank than at 
the stagnation point of this vehicle, so the volume of radiating gas is greatly increased. The 
vibrational temperature in this nonequilibrium flowheld is still high in this region, so the 
radiative emission is strong. 

The other variable of interest for coupled howhelds is the divergence of the radiative 
flux in the shock layer. To compare the two transport methods, Figs. 8 and 9 show contour 
plots of V ' <?R i n the computational domain (see Fig. 4), since it is difficult to see details in 
the physical domain. For a viscous how, this has the effect of stretching the boundary layer 
region, which extends to (/(max of 0.3 to 0.4 in these plots. The position of the first 0.0 
contour level in the boundary layer region is affected by different gridding in the two cases, 
and is not indicative of a difference in the solutions. The tangent slab result shown in Fig. 8 
is in fact just the derivative of ^r along each normal grid line, so the values at neighboring 
grid lines do not influence each other. The localized hot spot near £/£ ma x = 0.8 may be an 
artifact of this 1-D assumption. The grid problems in the shoulder area are evident from the 
distorted contours near £/£ ma x = 1, and in fact the last grid line is not included in the MDA 
solution. In all other aspects, the solutions from the two transport methods are qualitatively 
and often even quantitatively similar, giving confidence in the MDA method. 

The tangent slab calculation for this case required 7 CPU minutes on an IRIS 4D/340 
workstation. The MDA result took 35.8 CPU hours, or about 300 times longer. While this 
is a large penalty for the slight difference obtained in the results, several comments can be 
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made. Since the IRIS workstation has four processors, the MDA result required only about 
10 hours of real elapsed time. The speedup factor that could be obtained by running the 
method on the Cray is expected to be especially large, because it will allow the use of single 
rather than double precision arithmetic. Finally, the MDA algorithm is in its early stages, 
and significant reductions in computing time will probably be obtained as the algorithm 
matures. In particular, this result was obtained using a larger radiation subgrid (32 grid 
points normal to the body) than was suggested in Ref. 1 in order to obtain the best results 
for the evaluation of the method. Taking these factors into account, the MDA method has 
potential for providing accurate heating distributions in important cases, though it may 
never be used as commonly as the tangent slab method. More important, though, it can be 
used in flow regions (such as wakes) where the tangent slab method cannot be applied. 

Geosynchronous Return 

Results were also obtained at a flight condition representative of return from a geosyn- 
chronous orbit. The aerobrake configuration is a raked cone with a blunted elliptic nose 
(AFE). 19 Several studies of nonequilibrium radiative heating for this configuration have been 
made. 20-23 An axisymmetric LAURA flowheld solution in which this geometry was modeled 
by a sphere with an equivalent nose radius of 2.16 m 24 was generated for this study. 25 Re- 
sults will only be given for the stagnation region, since the sphere quickly deviates from the 
AFE geometry. The flowheld models used to obtain this solution are the baseline models 
of Ref. 26, including the use of the Park-87 chemical kinetics. The wall temperature was 
assumed to be constant at 1750 K. 

Radiation predictions were obtained for this case, again using both tangent slab and 
MDA transport methods. The tangent slab calculation, including calculation of the radiation 
properties k' v and j®, was completed in about 45 minutes of actual elapsed time on a Silicon 
Graphics 4D/320, while the MDA result required about 4 hours of elapsed time. Including 
the calculation of the radiation properties in the comparison means that the MDA method 
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requires just over 5 times more CPU time than the tangent slab method for this case. It is 
anticipated that the MDA run time can be further reduced as the code matures, as discussed 
above. 

The radiative heating rates predicted for this configuration at 9.3 km/s and 76 km alti- 
tude are shown in Fig. 10. These predictions are in the same range as earlier computations 9,27 
(note that different nose radii have been used for these computations). The MDA prediction 
is again lower than the tangent slab result, with the dip at the outflow boundary to be 
attributed to the outflow boundary condition used, rather than any physical effect. 

The radiative flux divergence for this flowheld is presented in Figs. 11 and 12 for the 
tangent slab and MDA transport, respectively. This time the physical domain is used, since 
details in the nose region are apparent in these coordinates. The approximations inherent 
in the tangent slab solution can be observed in this case as significant differences in the 
profiles along adjacent surface normal lines. As mentioned earlier, the tangent slab method 
employs a numerical differentiation to obtain dq^/dn while the MDA method requires no 
such differentiation. This may add to the error already incurred by ignoring the variation 
of the radiation properties in the direction parallel to the surface. This result suggests that 
for some flowheld cases, the MDA result might provide a smoother variation of the radiative 
flux divergence when howheld coupling is included, and thus might enhance the stability of 
coupled solutions. 

The amount of heating to geosynchronous return vehicles resulting from radiation in 
the ultraviolet (UV) portion of the spectrum has been the subject of some debate. The 
distribution of radiation predicted by the LORAN method is shown in Fig. 13, using both 
the tangent slab and MDA transport models. The relative importance of different spectral 
regions varies only slightly between the two transport methods. Both curves indicate that 
just over half of the radiative heating experienced by this vehicle results from the UV portion 
of the spectrum. This is a significant fraction of the total. Radiation in the UV spectral 
range is highly self-absorbed. If the amount of self-absorption is mispredicted only slightly, 
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the wall radiative flux may be significantly increased. 


Conclusions 

A Modified Differential Approximation (MDA) method for computing 3-D radiative 
transport in axisymmetric thermochemical nonequilibrium flows is developed. It employs the 
P-1 approximation to reduce the integro-differential governing equation of radiative transport 
to a set of partial differential equations. The numerical solution of these equations in a 
finite volume algorithm with real radiative property variations is discussed. The method is 
assessed in forebody flowhelds by comparison to the commonly used tangent slab transport 
approximation. The MDA method is intended for eventual use in flow areas, such as wakes, 
where the tangent slab approximation is known to be invalid. 

Predictions of the radiative heat flux to the forebody wall of an axisymmetric vehicle 
are obtained from the two transport methods for two representative flowhelds. Comparison 
of the results shows qualitative agreement, with wall radiative heating predictions from the 
MDA method about 10 to 20 percent lower than from the tangent slab. The differences 
increase away from the stagnation region due to geometric and howheld gradient effects 
not accounted for in the tangent slab approximation. The variation of the divergence of 
the radiative heat flux, which appears in the energy equation to couple the radiation and 
howheld properties, is also examined and found in some cases to be smoother in the MDA 
method. In addition to providing more accurate coupled results, this feature of the MDA 
method has the potential to enhance the stability of coupled solutions. 
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Figure Captions 


1. Typical High Speed Aerobrake Flowfield 

2. Flow Geometry and Radiation Boundary Conditions 

3. Variation of ~K V and j v for a Mars Return Flowfield 

4. Physical to Computational Domain Transformation 

5. Typical Convergence History for MDA Solution 

6. Wall Radiative Flux for Mars Return Case 

7. Comparison of MDA and Tangent Slab Result 

8. Radiative Flux Divergence in Computational Domain for Mars Return Case with 
Tangent Slab Transport 

9. Radiative Flux Divergence in Computational Domain for Mars Return Case with 
MDA Transport 

10. Wall Radiative Flux for Geosynchronous Return 

11. Radiative Flux Divergence for Geosynchronous Return with Tangent Slab Transport 

12. Radiative Flux Divergence for Geosynchronous Return with MDA Transport 

13. Spectral Distribution of Geosynchronous Return Wall Radiative Flux 
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Figure 1: Typical High Speed Aerobrake Flowfield 
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Figure 2: Flow Geometry and Radiation Boundary Conditions 
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Figure 3: Variation of and j v for a Mars Return Flowfield 
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Figure 6: Wall Radiative Flux for Mars Return Case 
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Figure 7: Comparison of MDA and Tangent Slab Result 
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Figure 8: Radiative Flux Divergence in Computational Domain for Mars Return Case with 
Tangent Slab Transport 
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Figure 9: Radiative Flux Divergence in Computational Domain for Mars Return Case with 
MDA Transport 
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Figure 10: Wall Radiative Flux for Geosynchronous Return 
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Figure 11: Radiative Flux Divergence for Geosynchronous Return with Tangent Slab Trans- 
port 
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Figure 12: Radiative Flux Divergence for Geosynchronous Return with MDA Transport 
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